Cite this Article
Qin Zhang-Rong, Chen Yan-Yan, Ling Feng-Ru, Meng Ling-Juan, Zhang Chao-Ying. A mass-conserved multiphase lattice Boltzmann method based on high-order difference. Chinese Physics B, 2020, 29(3): 034701  
Permissions
A mass-conserved multiphase lattice Boltzmann method based on high-order difference
Qin Zhang-Rong, Chen Yan-Yan, Ling Feng-Ru, Meng Ling-Juan, Zhang Chao-Ying †
Guangxi Collaborative Innovation Center of Multi-source Information Integration and Intelligent Processing, Guangxi Normal University, Guilin 541004, China

 

† Corresponding author. E-mail: zhangcy@gxnu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11862003 and 81860635), the Key Project of the Natural Science Foundation of Guangxi Zhuang Autonomous Region, China (Grant No. 2017GXNSFDA198038), the Project of Natural Science Foundation of Guangxi Zhuang Autonomous Region, China (Grant No. 2018GXNSFAA281302), the Project for Promotion of Young and Middle-aged Teachers’ Basic Scientific Research Ability in Guangxi Universities, China (Grant No. 2019KY0084), and the “Bagui Scholar” Teams for Innovation and Research Project of Guangxi Zhuang Autonomous Region, China, and the Graduate Innovation Program of Guangxi Normal University, China (Grant No. JXYJSKT-2019-007).

Abstract

The Z–S–C multiphase lattice Boltzmann model [Zheng, Shu, and Chew (ZSC), J. Comput. Phys. 218, 353 (2006)] is favored due to its good stability, high efficiency, and large density ratio. However, in terms of mass conservation, this model is not satisfactory during the simulation computations. In this paper, a mass correction is introduced into the ZSC model to make up the mass leakage, while a high-order difference is used to calculate the gradient of the order parameter to improve the accuracy. To verify the improved model, several three-dimensional multiphase flow simulations are carried out, including a bubble in a stationary flow, the merging of two bubbles, and the bubble rising under buoyancy. The numerical simulations show that the results from the present model are in good agreement with those from previous experiments and simulations. The present model not only retains the good properties of the original ZSC model, but also achieves the mass conservation and higher accuracy.

PACS: ;47.11.-j;;05.50.+q;;47.61.Jd;;02.70.Bf;
Keyword:;lattice Boltzmann method;;high-order difference;;mass conservation;;large density ratio;
1. Introduction

The multiphase flow problems are ubiquitous in nature, industrial and agricultural productions. It is of great significance to study the basic theory and law of multiphase flow for people’s production and living and promoting the development of industry and agriculture. In order to study these problems numerically, lattice Boltzmann method (LBM) has been developed into a powerful tool. Due to numerous advantages of the LBM, including solid physical foundations, high efficiency, complete parallelization, and easy processing of complex boundaries, it has been successfully used in the numerical simulation of complex multiphase fluid systems.[1–3]

In the past decades, many efforts have been made to develop multiphase lattice Boltzmann (LB) models, and several main popular models are proposed, i.e., the color-gradient model proposed by Gunstensen et al.,[4] the pseudo-potential model (SC model) by Shan and Chen,[5,6] the free energy model by Swift et al.,[7] and the HCZ model by He et al.[8] Among them, the free energy model has been favored by many researchers because it can solve the problem of interface tracing, restore the Cahn–Hilliard (CH) equation, satisfy the local mass and momentum conservation, conform with the thermodynamic theory, etc. However, the numerical instability will be caused by the original free energy model with large density ratio, and the Galilean invariance is not satisfied when there is a large density gradient near the interface.[9] To solve this problem, Inamuro et al.[10] proposed a projection method for the free-energy model to achieve a large density ratio, in which the Poisson equation for pressure convergence was used. However, solving the Poisson equation would spoil the simplicity and efficiency of the LBM. Lee and Lin[11] later developed a stable discretization scheme (LL model) for two-phase flows with high density and viscosity ratio, in which the pressure gradient term is discretized by different ways before and after streaming step, which makes the scheme implementation relatively complicated and the amount of computation large. Zheng et al.[12] presented a multiphase LB model with large density ratio (ZSC model) based on the free-energy model, in which the CH equation is accurately recovered without any additional terms and the Galilean invariance property remains unchanged. In addition, owing to the small change of the average density in this model, it is very stable and efficient, and more importantly, it can be employed to simulate the multiphase flows with large density ratio beyond 1000. Due to good performance of the ZSC model, it has been favored by many researchers. Huang et al.[13] extended the ZSC model to three-dimensions to simulate a bubble rising phenomenon. Cheng et al.[14] simulated the multiple bubbles rising under buoyancy in a quiescent viscous incompressible fluid with a three-dimensional (3D) ZSC model. Recently, He et al.[15] employed the ZSC model to simulate and analyze the role of wettability and surface tension in forming the ink-jet printing droplets. However, owing to the presence of the diffusion term and the numerical dissipation caused by the discretization of convection term in CH equation, the mass of each phase in the ZSC model cannot be conserved exactly.[9] van der Sman and van der Graaf et al.[16] found that the droplet mass under shear decreases continuously in the evolution process. Huang et al.[13] showed that the volume of a 3D bubble rising under buoyancy dwindles over time when Reynolds numbers is high, and in an extreme case, the bubble volume loss can increase up to 75%. Zheng et al.[17] also illustrated that a non-physical shrinkage of the bubble may be produced in the diffuse interface method. In order to solve the non-conservative mass problem in phase-interface diffusion method, Son[18] introduced a volume correction equation for the order parameter in the level set formulation to keep the mass conservation, which needs multiple iterations to correct the order parameter in each time step. Based on this method and considering the effect of density changes, Chao et al.[19] recasted the volume correction equation and presented a filter-based, mass-conserved LB model for improving the HCZ model, in which high density ratio (up to 100) and mass conservation are achieved. By using Chao et al.’s mass correction method and an effective surface tension formula, Huang et al.[20] developed a mass-conserved axisymmetric multiphase LB model of simulating the bubble rising, which retains mass conservation of the bubble, and their the maximum density ratio is 15.5. Later, Wang et al.[9] considered that the mass loss or increase in a phase state goes through the interface zone by diffusion, and they proposed a mass correction method to compensate for the mass loss or offset the mass increase in the volume of diffuse layer to keep the mass conservation at each time step. Their results showed that the mass is well-conserved with high density ratio. Recently, Niu et al.[21] also presented a multiphase LB model to improve Shao et al.’s model[22] by using a mass correction method similar to that in Wang et al.’s work.[9] The mass conservation for each phase is maitained, but their density ratio is quite small. In brief, the above two correction methods are both favored by researchers. The volume correction method is relatively simple, but the efficiency is affected by multiple iterations in each correction. The mass correction method based on the phase interface is more efficient, but it is not suitable for the ZSC model which cannot accurately obtain the local density by the order parameter.

On the other hand, the solutions to the gradient of the order parameter are involved in the calculations of chemical potential and velocity in the ZSC model, and the accuracy for calculating these gradients is particularly important. The central difference method is often used to calculate the numerical gradient in multiphase models.[11,14,23] In general, this method is sufficiently accurate. However, in the diffusion interface method with high density ratio, the results calculated by central difference method may deviate greatly from the theoretical solutions, which may lead the numerical simulation to be unstable. Therefore, the introduction of a high accurate difference method instead of traditional central difference method is expected to be able to further improve the computational accuracy of the ZSC model.

Through the above investigations, we know that although great progresses have been made in the numerical study of multiphase flow based on free energy model, solving the problems of large density ratio and mass non-conservation in the model are still an important research subject. In this paper, we are to improve the original ZSC model in terms of the mass conservation and computational accuracy in order to develop a mass-conserved and high accuracy multiphase LB model for simulating multiphase flows, in which the high-order difference is introduced to calculate the gradient of the order parameter to improve accuracy and a mass correction method proposed by Chao et al.[19] is utilized to solve the non-conservative mass problem. As a result, the present model not only enjoys the advantages of the original model such as large density ratio, stability, efficiency and Galilean invariance, but also can solve its non-conservative mass problem and further improve its accuracy. We will investigate the numerical methods of calculating the gradient of the order parameter and reveal the importance of its numerical accuracy, and demonstrate the performances of the improved model by testing several cases such as a bubble in a stationary flow, the merging of two bubbles, and the bubble rising under buoyancy.

The rest of this paper is organized as follows. Firstly, a brief introduction is given to the original ZSC model. Secondly, the central difference and the high-order difference methods are investigated. Thirdly, a mass-conserved multiphase LB model based on the original model is constructed by introducing a mass correction term. Fourthly, the performances of the present model are investigated through using several testing cases. Finally, a brief conclusion is drawn.

2. Method
2.1. Macroscopic governing equations

In this study, two-phase fluid with different densities is considered. The high density and low density are represented by ρH and ρL respectively. The flow can be described by the Navier–Stokes (NS) equations and a Cahn–Hilliard equation as[24–26]

where p is the pressure tensor, θM is the mobility, Fb is the body force, μϕ is the chemical potential, n and ϕ are defined as n = (ρA + ρB)/2, ϕ = (ρA - ρB)/2,[26] respectively, ϕ is the order parameter indicating the two phases, while ρA and ρB are the densities of fluid A and fluid B, respectively. Depending on the initial conditions, they may be ρL or ρH.

In Eq. (2), the term ∇ ⋅ p relates to the surface tension, which can be written as a potential term[25,27]

where and cs is the sound speed.

In combination with Eq. (4), the momentum equation (2) can be rewritten as

The chemical potential μϕ can be derived from the free energy density function, and a free energy function in a closed volume with a mixture of two fluids is adopted as[25,28,29]

where V is a control volume, ψ is the bulk free energy density per unit mass for the homogeneous system, κ is a coefficient relating to the surface tension and the interface layer thickness. The square of gradient term relates to the variation of the density and contributes to the excess of the free energy in the interface, which determines the surface energy.[28] It is chosen as a double-well form

where A is an amplitude parameter, which is used to control the interaction energy between the two phases. This form will contribute to two equilibrium states, ϕ* and −ϕ*. The chemical potential μϕ is calculated from[13]

where A = 3 σ(4wϕ*4), κ = A(wϕ*)2/2,[13] σ is the surface tension coefficient, W is the interface thickness, and the expected order parameter is ϕ* = (pH – pL)/2.

2.2. Interface capturing LB equation

The interface capturing is modeled by the convective CH equation (3). The distribution function gi is used to track the interface between two phases in each lattice node. To solve Eq. (3), a modified LB equation is adopted[30]

where the collision operator is , with gi being the distribution function for tracking interface, and τϕ being the dimensionless relaxation time, ci is the lattice velocity, and here the 3D D3Q7 velocity model is used, q = 1/(τϕ + 0.5) is a constant coefficient, is the equilibrium distribution function, and it is in the following form:

the coefficients can be chosen as

where Γ is the mobility. The macroscopic variable (order parameter) ϕ is calculated from

2.3. Continuity and momentum LB equation

The isothermal incompressible Newtonian fluid is modeled by Eqs. (1) and (2). Here, another distribution function fi is utilized for simulating the fluid flow in each lattice node, and in order to solve Eqs. (1) and (2), the second LB equation is employed as[12]

where fi is the distribution function for fluid flow, the lattice velocity ci takes the 3D D3Q19 velocity model, Si is a source term added into LB equation to mimic the body force term of the NS equations, and it is given as

the equilibrium distribution function is taken as

where the coefficients are chosen as[13]

The macroscopic variable n and velocity u are calculated from

2.4. Central difference and high-order difference method

It can be seen from Eqs. (8) and (16) that the second-order and the first-order gradient of the order parameter ϕ are involved in solving the chemical potential μϕ and velocity u respectively. The performance of numerical method of calculating gradient is likely to affect the accuracy and stability of the multiphase flow model. In this subsection, we will investigate the accuracy of the central difference method (CDM) and high-order difference method (HDM) in some specific cases, and then propose a better numerical scheme to improve the accuracy of the ZSC model.

2.4.1. Central difference method

Consider a continuous and derivable function y = f(x) [a,b],a ≤ b, and the interval [a,b] is equidistantly divided into n sub-intervals with the spacing h = (b – a)/n and separated by xi(i = 0,1,n – 1,n). According to the principle of CDM, the numerical derivative at xi can be solved by

From Eq. (17), it can be seen that the derivate at node xi is determined by its neighbouring function value of f(xi – 1) and f(xi + 1), and independent of f(xi) at its own position. In the actual simulations of multiphase flows, a weighted central difference scheme is usually adopted to calculate the gradient. For the D3Q19 model, the numerical gradient of the order parameter can be evaluated from[11]

For a phase transition system, its equilibrium distribution of density or order parameter can be usually approximated into a hyperbolic tangent function. Here, consider a bubble in the center of the stationary flow field, then the hyperbolic tangent function will be used to initialize the flow field and written as follows:

where the computational domain is DX × DY × DZ = 110 × 110 × 110 in lattice unit. The initial position of the bubble is (x0,y0,z0) = (DX/2,DY/2,DZ/2), and its radius is R = 10. The liquid density takes ρH = 1000, the gas density takes ρL = 1, and the density ratio is 1000. The initial order parameter is ϕ* = (ρH – ρL)/2,–ϕ* for the bubble and ϕ* elsewhere. The surface tension is σ = 0.1, the interface thickness is W = 5, the relaxation times are τn = 0.7 and τϕ = 0.7 respectively, the migration coefficient is Γ = 100, and the periodic boundary conditions are applied to all boundaries. Unless otherwise specified, all the parameters in other testing cases in this paper are the same as those above.

The first and second derivatives of the order parameter calculated by CDM across the interface of the bubble are shown in Fig. 1, the corresponding exact results are also plotted for comparison, and their absolute errors are shown in Fig. 2. In order to investigate the errors from CDM with different interface thickness, the standard errors between the derivatives calculated by CDM and the exact results are shown in Fig. 3. From Figs. 1 and 2, it can be clearly seen that there exists a significant discrepancy between the numerical result calculated by CDM and the theoretical value, and the discrepancies of the second derivative are greater than those of the first one. Figure 3 shows that the deviations increase sharply with interface thickness W decreasing. These discrepancies can affect the simulation accuracy or even produce incorrect results in the numerical simulation of multiphase flow.

Fig. 1. Computational results from CDM and HDM, showing (a) the first derivative (in x-axis direction) and (b) the second derivative.
Fig. 2. Absolute errors from CDM and HDM for (a) the first derivative (in x-axis direction) and (b) the second derivative.
Fig. 3. Plots of standard error versus interface thickness from CDM and HDM for (a) the first derivative (in x-axis direction) and (b) the second derivative.
2.4.2. High-order difference method

Considering that the properties of the fluid nodes near the node to be solved are similar, the derivation using more information about adjacent nodes should obtain more accurate results. We construct a multi-order difference method by performing a Taylor series expansion at different neighbor nodes to build polynomials, and then evaluate the correlation coefficients of the polynomials with a method of undetermined coefficients. In order to give consideration to the accuracy and efficiency of this method, a five-point higher-order difference method is used in this study[31,32]

To verify the performance of HDM, the similar calculations in Subsubsection 2.4.1 with HDM are also performed. The corresponding results are also plotted in Figs. 1–3, respectively. As shown in Fig. 1, the derivatives from HDM are closer to the exact values than the ones from CDM. Figure 2 obviously shows that the accuracy of HDM is much higher than that of CDM according to the absolute error of the derivative. Figure 3 shows that as the interface thickness decreases, the standard error from HDM increases, but the error level is much lower than that from CDM. In order to investigate the computational efficiency of CDM and HDM, a case of a bubble in a stationary flow is tested. The bubble radius is R = 30, and the computational domains are 100 × 100 × 100, 150 × 150 × 150, 200 × 200 × 200, and 300 × 300 × 300, respectively. The CDM and HDM are used to calculate separately the first and second derivatives of the order parameter in Eqs. (8) and (16). The computing platform is Windows 10 64 bit, Intel(R) Core(TM) i9-9900K CPU, and 32 GB RAM. The time spent in running the program with 10000 steps is shown in Table 1. As indicated in Table 1, the time spent in using HDM is much less than that in using the CDM, and the average computational efficiency of HDM is about 1.53 times higher than that of CDM. The main reason for the above result is that the using of HDM requires relatively less lattice nodes. For the first derivative of the order parameter in each direction, the using of HDM (see Eq. (23)) only needs to consider 4 neighbor nodes in the axial direction, while the using of CDM (see Eqs. (18)–(20) needs to consider the adjacent 10 neighbor nodes. For the second-order derivative of the order parameter, the using of HDM (see Eq. (24)) needs to employ 15 nodes in the x-, y-, and z-axis directions. However, the using of CDM (see Eq. (21)) needs to use 18 neighbor nodes. Through the above analysis, we can conclude that derivatives calculated with HDM may be more accurate and efficient than the ones calculated with CDM in simulations. Therefore, the five-point high-order difference method will be introduced to calculate the gradient of the order parameter in Eqs. (8) and (16).

Table 1

Comparison of computational efficiency between CDM and HDM.

.
2.5. Mass conservation correction

Mass conservation is an important problem in the numerical simulation of multiphase flow. Although the discretization of the CH equation has the form of mass conservation, Ding et al.[33] pointed out that the mass conservation does not mean that the volume enclosed by any given contour remains unchanged. In order to solve the problem of non-conservative mass in the original ZSC model, a volume conservation correction method proposed by Chao et al.[19] is adopted

where V is the volume of the gas or liquid phase before the correction after streaming step, V0 is its initial volume. After each streaming step, at every artificial time τ, we calculate the volume difference Vd = |V–V0|, if Vd > err (an error criterion defined by user), update the order parameter ϕ(x, t + Δt) with the conservation corrected values calculated from Eqs. (25) and (26), and then recalculate the Vd according to the updated order parameter. If Vd is still more than the err, we iterate the above processes until Vd is less than the err or the iteration steps reach the maximum iteration step, and then continue and enter into the next LB evolution iteration step. In order to improve the computational accuracy of Eqs. (25) and (26), the five-point higher-order difference method is used. The artificial time step ∂τ = 0.15/V0 is adopted in our simulation.[20]

3. Numerical verification
3.1. Bubble in stationary flow
3.1.1. Laplace law

The Laplace is a basic testing case to verify the surface tension property of a multiphase flow system. According to this law, when a bubble and its adjacent liquid reach a stable state, the relationship between the pressure difference across the bubble interface and the surface tension is as follows:[34,35]

where Pin and Pout are the pressure inside and outside the bubble respectively, k = 2 is for the 3D model, σ is the surface tension, and R is the radius of the bubble. A spherical bubble is placed in the center of a flow field with a mesh size of 130 × 130 × 130. The surface tension is σ = 0.1, and the interface thickness is W = 4. The radius of the bubble is taken as 15, 20, 25, 30, 35, 40, 45, and 50, separately. The period boundary condition is imposed on the all sides.

Figure 4(a) shows the relationship between the pressure difference and the radius of the bubble, and it can be seen that the results from the HDM are in excellent agreement with those from the Laplace law. With the decrease of the radius, the maximum absolute error between the results from CDM and theoretical solutions reaches up to 10%. In addition, figure 4(b) shows the effect of the interface thickness on the surface tension. It can be easily observed that the surface tension obtained by the CDM deviates significantly from the theoretical value as the thickness of the interface decreases when the thickness of the interface is smaller than 4.5. However, the surface tension obtained by the HDM well matches the theoretical value when the thickness of the interface is bigger than 2.5, and it just decreases a little when the thickness of the interface is smaller than 2.5. Form Figs. 4(a) and 4(b), it can also be seen that evenintroducing mass correction into the CDM cannot improve the computational accuracy.

Fig. 4. (a) Pressure differences versus 1/R and (b) surface tension versus interfacial thickness (with radius of bubble R = 20).
3.1.2. Mass conservation tests

In order to investigate the mass conservation of the improved model, the mass conservation are tested by using the CDM and HDM with or without mass correction in a 3D stationary flow. Figure 5 shows the evolution of the bubble mass with radius in different times. It can be seen that the bubble mass obtained by the original ZSC model with the CDM or the HDM decreases continuously with the increase of the number of evolution steps. When the number of evolution steps is 50000, the radius of bubble is 20, 30, 40, and 50, the mass loss is 1.726%, 0.638%, 0.368%, and 0.239% for CDM respectively, and the corresponding mass losses for the HDM are almost the same. However, the bubble mass obtained by the CDM and the HDM with mass correction are both maintained very well during the whole simulation.

Fig. 5. Evolution of bubble mass with time for different radii.
3.2. Two merging bubbles

Bubbles merging is a typical phenomenon of bubble motion. In order to study the performance of the present model with mass correction under the motion of bubbles, a testing case of two merging bubbles is executed. In this simulation, the surface tension is σ = 1.0 and the initial spacing between two bubbles is 5. Figure 6 shows the time evolution of two bubbles merging, and the simulation results and processes are consistent with those in Refs. [9] and [12]. Figure 7 displays the mass changes of the two merging bubbles during the evolution. It can be seen that if there is no mass correction, with the two bubbles merging slowly, the mass of the two merging bubbles gradually decreases with the evolution step increasing, and when the number of evolution steps is 200000, the mass loss of the bubbles is 6.75% for CDM and 6.69% for HDM. By comparing the results in Subsubsection 3.1.2, we also find that the mass loss of dynamic bubble is greater than that of static bubble. However, the results indicate that the mass losses of the two merging bubbles are almost zero when the mass correction is introduced into the original model with the CDM or HDM.

Fig. 6. Evolutions of two bubbles merging with time.
Fig. 7. Evolutions of mass with time for two merging bubbles.
3.3. Single bubble rising

The complex motion characteristics such as deformation and oscillation will occur during a bubble rises under buoyancy, the deformation of the bubble relates to surface tension, and the deformation degree relates to Reynolds (Re) number, Eotvos (Eo) number, and Morton (Mo) number, and they can be calculated from the following formulas:[36,37]

The dimensionless parameter for the evolution time is calculated from[14]

where t is evolution step, and the bubble velocity can be calculated from[12]

where U is the bubble velocity component in the vertical direction.

In order to further verify the performance of the present model (HDM with mass correction) under the complex motion conditions, a case of single bubble rising will be teseted. Initially, a spherical bubble is located at the bottom of the simulation domain with 120 × 120 × 240, 30 lattice units high away from the bottom wall, and its diameter is d = 20. The density ratio between liquid and gas is 1000 and viscosity ratio is 100, the interface thickness is W = 4, and the migration coefficient is Γ = 1000. The extrapolation boundary condition is applied to the top boundary, and the half bounce-back boundary conditions are imposed on the other boundaries.

Table 2 gives the comparison between the Re numbers from the present model, CDM with mass correction and other results from Huang et al.[13] and the experiments of Clift et al.[36,37] in different Eo and Mo numbers. It can be seen that the results from Huang et al. are quite close to the experimental results from Clift et al. in the cases C2 and C3, except for the case C5, and its relative error is up to 16.853%. The maximum relative error from the present model is 8.964% in the case C4, and the minimum one from the present model is only 0.172% in the case C2. Meanwhile, it can also be seen that the relative error of the present model is smaller than that of the CDM with mass correction. Therefore, the present model has higher accuracy, and its results are in good agreement with the experimental results in all the cases.

Figure 8 illustrates the evolutions of Re number with the evolution time in some cases (C1–C4) in detail. It can be seen that with the increase of the evolution time, the Re number increases rapidly then decreases slowly, and finally reaches a steady state for each of all the cases. Moreover, the larger the Re number, the longer the time it takes to reach a steady state and the greater the discrepancy of Re number between the experiment and the present model. We also find that the Re number from HDM with mass correction is better than those from other methods, andit is closer to the experimental result for each of all the cases.

Fig. 8. Plots of evolution of Re number with time in four cases (C1–C4).
Table 2

Comparison between Re numbers from the present model and other results with different Eo and Mo numbers.

.

Figure 9 shows the evolution of the rising bubble mass in four cases (C1)–(C4) in the Table 2. It is obvious that without the mass correction, the rising bubble mass from CDM and HDM models decreases slowly with the evolution time increasing for each of all the cases, and the larger the Re number, the larger the mass loss is. Nevertheless, the mass loss from the HDM is much smaller than that from the CDM. If the HDM or CDM with mass correction is applied to the original model, the rising bubble mass is conserved very well in all evolution time. By analyzing the above testing cases, including the bubble in a stationary flow, the merging of two bubbles, and the bubble rising under buoyancy, we can draw a conclusion that the mass conservation can be achieved by introducing the mass correction method into the original model, and the main contribution of mass conservation is the mass correction method, which is unrelated to whether CDM or HDM is used in the model. However, by using the HDM the computational accuracy and efficiency can be improved.

Fig. 9. Evolutions of rising bubble mass with time in four cases (C1)–(C4).

Besides the comparison of the Re number and the mass conservation, to verify the performance of the present model, the bubble shapes under different parameters are also studied. The bubble shape is affected by surface tension, and Eo number is related to surface tension. A high Eo number means a low surface tension, which can make bubble greatly deformed. In addition, Re number relates to the deformation in the bubble’s vertical direction, and the larger the Re number, the larger the deformation will be. In general, the final bubble shape is determined by both Re and Eo numbers.[20] Figure 10 shows the comparison between bubble shapes from the present model and those from Huang et al.[13] and experiment.[36,37] For case C1, due to small Eo and Re, the final bubble shape remains almost unchanged when it rises. As Eo or Re increases, the final bubble shape becomes an oblate ellipsoidal cap in case C3, an oblate ellipsoidal (disk) in case C5, and a spherical cap (closed wake) in case C6 respectively. These shapes from the present model all consistent well with the experimental ones.[36,37] The final bubble shapes from Huang et al.[13] are also very close to those from the experiments.[36,37] in cases C1, C3 and C5. However, they found that the bubble deforms and rises in an oscillatory manner, which has some deviation from the experimental results.[13] In order to further study the effect of Mo number on bubble shape, a comparison among bubble shapes from different methods is made with Eo = 116 and different Mo numbers. As shown in Fig. 11, with the Mo number decreasing, their differences in deformation at the bottom of the spherical crown increase gradually. The simulation results from the present model are consistent with those from Hua et al.,[38] Cheng et al.,[14] Huang et al.,[13] Ren et al.,[39] Li et al.,[40] and the experiment.[36]

Fig. 10. Comparison among bubble shapes from the present model, Huang et al.,[13] and experiments.[36,37]
Fig. 11. Comparison among bubble shapes with Eo = 116 and different Mo numbers by using different methods.
4. Conclusions

In this paper, a mass-conserved and high accuracy multiphase LB model based on the ZSC model is developed, in which the high-order difference is introduced to calculate the gradient of the order parameter and thus improving the accuracy, and a mass correction method is employed to solve the non-conservative mass problem of the original model. In order to verify the performance of the present model, several cases such as a bubble in a stationary flow, the merging of two bubbles, and the bubble rising under buoyancy are tested. Some conclusions are drawn from the present study as follows.

(i) Whether we use simple single static bubble simulation or bubble rising simulation with complex deformation motion characteristics, in all testing cases, through introducing the mass correction method into the original model, the bubble mass is conserved very well, and the problem about non-conservation of mass in the original modelis solved.

(ii) By introducing HDM into the original model, the simulation results from HDM are much better than those from CDM in all the testing cases, the accuracy of the present model is improved.

(iii) In the simulations of bubble rising, the Re numbers and the final bubble shapes from the present model are compared with those from the experiments and other models or methods. The results from the present model are in good agreement with the experimental ones, and better than those from other models and methods mentioned when the Re number is relatively small.

(iv) In all our simulations, the density ratio for each of two phases is 1000, the simulations are very stable, and the present model keeps the advantages of the original ZSC model, such as large density ratio and stability very well.

Due to its high accuracy, mass conservation and large density ratio, the present multiphase model can be expected to be applied reliably to more general complex fluid systems and obtain some good results. Nevertheless, in the case of bubble rising with high Re number, the instability for the velocity and bubble shape may appear, and our near future work is to further improve the present model to solve this problem.

Reference
[1] Ye H F Kuang B Yang Y H 2019 Chin. Phys. 28 014701
[2] Wen B H Chen Y Y Zhang R L Zhang C Y Fang H P 2013 Chin. Phy. Lett. 30 064701
[3] Hu M D Zhang Q Y Sun D K Zhu M F 2019 Acta Phys. Sin. 68 030501 in Chinese
[4] Gunstensen A K Rothman D H Zaleski S Zanetti G 1991 Phys. Rev. 43 4320
[5] Shan X Chen H 1993 Phys. Rev. 47 1815
[6] Shan X Chen H 1994 Phys. Rev. 49 2941
[7] Swift M R Osborn W R Yeomans J M 1995 Phys. Rev. Lett. 75 830
[8] He X Y Chen S Y Zhang R Y 1999 J. Comput. Phys. 152 642
[9] Wang Y Shu C Shao J Y Wu J Niu X D 2015 J. Comput. Phys. 290 336
[10] Inamuro T Konishi N Ogino F 2000 Comput. Phys. Commun. 129 32
[11] Lee T Lin C L 2005 J. Comput. Phys. 206 16
[12] Zheng H W Shu C Chew Y T 2006 J. Comput. Phys. 218 353
[13] Huang H Zheng H Lu X Shu C 2010 Int. J. Numer. Meth. Fluids 63 119
[14] Cheng M Hua J Lou J 2010 Comput. Fluids 39 260
[15] He B Yang S C Qin Z R Wen B H Zhang C Y 2017 Sci. Rep. 7 11841
[16] Van Der Sman R G M Van Der Graaf S 2008 Comput. Phys. Commun. 178 492
[17] Zheng L Lee T Guo Z David R 2014 Phys. Rev. 89 033302
[18] Son G 2001 Num. Heat Tr. B-Fund. 39 509
[19] Chao J Mei R Singh R Wei S 2011 J. Num. Meth. Fluid 66 622
[20] Huang H Huang J Lu X 2014 J. Comput. Phys. 269 386
[21] Niu X D Li Y Ma Y R 2018 Phys. Fluids 30 013302
[22] Shao J Y Shu C Huang H B Chew Y T 2014 Phys. Rev. 89 033309
[23] Kupershtokh A L Medvedev D A Karpov D I 2009 Comput. Math. Appl. 58 965
[24] Jacqmin D 1999 J. Comput. Phys. 155 96
[25] Kendon V M Cates M E Desplat J Pagonabarraga I Bladon P 2001 J. Fluid Mech. 440 147
[26] Takada N Misawa M Tomiyama A Hosokawa S 2001 J. Nucl. Sci. Technol. 38 330
[27] Jamet D Torres D Brackbill J U 2002 J. Comput. Phys. 182 262
[28] Rowlinson J S Widom B 1989 Molecular Theory of Capillarity Oxfordshire Clarendon Press 356
[29] Jamet D Lebaigue O Coutris N Delhaye J M 2001 J. Comput. Phys. 169 624
[30] Lamura A Succi S 2003 Int. J. Mod. Phys. 17 145
[31] Wang Y 2011 Mathematics in Practice and Theory 41 16 in Chinese
[32] Wang Y 2009 J. Tianjin Univ. Technol. 25 37 in Chinese https://www.cnki.com.cn/Article/CJFDTotal-TEAR200904010.htm
[33] Ding H Spelt P D M Shu C 2007 J. Comput. Phys. 226 2078
[34] Batchelor G K 1967 An introduction to fluid dynamics London Cambridge University Press 658
[35] Premnath K N Abraham J 2007 J. Comput. Phys. 224 539
[36] Bhaga D Weber M E 1981 J. Fluid Mech. 105 61
[37] Clift R Grace J R Weber M E 2005 Bubbles, drops, and particles New York Dover Publications 381
[38] Hua J Stene J F Ping L 2008 J. Comput. Phys. 227 3358
[39] Ren F Song B Sukop M C 2016 Ad. Water Resour. 97 100
[40] Li X Gao D Y Hou B L 2019 Chem. Eng. Sci. 193 76
[41] Hua J Lou J 2007 J. Comput. Phys. 222 769